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Jarzynski's identity for the free energy difference between two equilibrium states can be viewed 
as a special case of a more general procedure based on phase space mappings. Solving a system's 
equation of motion by approximate means generates a mapping that is perfectly valid for this pur- 
pose, regardless of how closely the solution mimics true time evolution. We exploit this fact, using 
crudely dynamical trajectories to compute free energy differences that are in principle exact. Nu- 
merical simulations show that Newton's equation can be discretized to low order over very large time 
steps (limited only by the computer's ability to represent resulting values of dynamical variables) 
without sacrificing thermodynamic accuracy. For computing the reversible work required to move a 
particle through a dense liquid, these calculations are more efficient than conventional fast switching 
simulations by more than an order of magnitude. We also explore consequences of the phase space 
mapping perspective for systems at equilibrium, deriving an exact expression for the statistics of 
energy fluctuations in simulated conservative systems. 



I. INTRODUCTION 

The maximum work theorem, a consequence of the sec- 
ond law of thermodynamics, states that the amount of 
work performed by a thermodynamic system during a 
transformation from a specific initial state A to a spe- 
cific final state B is less than energy difference between 
the two states 0. The work W is maximum and equal 
to the free energy difference, or reversible work, if the 
transformation is carried out reversibly. Equivalently, 
the average work performed on a system during such a 
transformation is bounded from below by the free energy 
difference AF, 



(W) > AF. 



(1) 



The notation (• • •) implies an average over many in gen- 
eral irreversible transformations initiated in an equilib- 
rium state. (For macroscopic systems every individual 
transformation will require the same amount of work but 
for small systems work fluctuations occur.) Remarkably, 
the inequality can be turned into an equality by con- 
sidering exponential averages 



exp(-/3Ai^) = (cxp(-/JVK)). 



(2) 



where (3 = l/k^T is the inverse temperature and fee is 
Boltzmann's constant. This identity, proven by Jarzynski 
and later by Crooks 01 under very general conditions, 
relates the statistics of irreversible work to equilibrium 
free energy differences. 

The Jarzynski identity can be used to calculate free 
energy differences in computer simulations of molecular 
systems 0, 0, IS • In statistical mechanical terms the 
free energy difference AF = Fb ~ Fa between a system 
at temperature T with Hamiltonian Ti.B{x) and another 
one at the same temperature with Hamiltonian Ti.A{x) is 



given by 



AF = -fcRTln 



/ dx exp{— /3Wb(x)} 
/ dx exp{— /37i^(x)} 



--fcBTln^ 



(3) 

where the integration extends over the entire phase space 
and X = {q, p} includes the positions q and momenta p of 
all particles. In the above equation, Qa and Qb are the 
canonical partitions functions of systems A and B. To 
calculate AF using Jarzynski's identity we introduce a 
parameter dependent Hamiltonian Ti.(x, A) defined such 
that Ha and Ti. b arc obtained for particular values of the 
control parameter, T-1(x,Xa) = Ti.Aix) and T-1{x,Xb) ~ 
T-Lb{x). By switching the control parameter A from Xa 
to Xb wc can continuously transform T-Ia{x) into TLb{x). 
If this is done over a time t, while the system evolves 
from particular initial conditions .tq , the work performed 
on the system is 



W = 



(4) 



Its value depends on the initial conditions xq and on 
the particular way the control parameter A(t) is switched 
from its initial to its final value. According to Jarzyn- 
ski's identity, the free energy difference can be evaluated 
by averaging the work exponential e~^^ over many such 
transformations. Specifically, this average is performed 
over a canonical distribution of initial conditions in the 
initial equilibrium state, 



exp(-/3AF) = j dxo p{xo) exp{- PW (xq)} , 



(5) 



where p{xo) = cxp{~ (3TiA{x)} / Q a- 

In fast switching simulations based on the Jarzynski 
identity, non-equilibrium trajectories are generated by 
approximately integrating the equation of motion, typ- 
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ically through a truncated Taylor expansion of the time- 
evolving phase space point x{t). The fidelity of trajecto- 
ries obtained in this way to true microscopic dynamics is 
determined by the time interval over which a low-order 
Taylor expansion is assumed to be accurate. Usually, 
the time step is chosen to be small, so that the total 
energy is nearly conserved when control parameters are 
held constant (in an isolated system) [a,[l3 • this paper 
we show that fast switching trajectories integrated with 
large time steps, while perhaps poor simulations of dy- 
namics, suffice to compute exact free energy differences. 
This new approach, which can increase the efficiency of 
fast switching simulations by up to two orders of magni- 
tude, is based on a generalization of Jarzynski's identity 
for general phase space mappings [Til ITa |. Jarzynski's 
original expression corresponds to the particular phase 
space mapping provided by the dynamical propagator. 
His result is valid, however, for any invcrtiblc phase space 
mapping. One could just as well use a concatenation 
of highly approximate molecular dynamics steps, the re- 
sult of integrating equations of motion to low order over 
large time intervals, to map points in phase space. Al- 
though such large time step trajectories are not accurate 
dynamical pathways, expressions for the free energy re- 
main exact. Due to the reduced cost of large time step 
trajectories, a considerable efficiency increase is possible. 

The remainder of the paper is organized as follows. 
The formalism and the justification of the large time step 
approach arc presented in section ^ The efficiency of 
the resulting algorithm is discussed in section IIIII In 
section Hvl we demonstrate the validity of this algorithm 
by calculating the reversible work to transform a simple 
one-dimensional energy landscape, and that to drag a 
particle through a Lennard-Jones fluid. Conclusions are 
given in section IVll 



II. FORMALISM 

A. Jarzynski's identity for phase space mappings 

The deterministic time evolution of a classical many- 
particle system can be viewed as mapping every point 
in phase space to another: a system initially at a;o will 
be located at Xt = (t>t{xo) after a time t. The function 
(j)t is called the propagator of the system. Since the sys- 
tem evolves deterministically the point Xt is completely 
determined by the initial conditions xq. The time re- 
versibility of equations of motion further ensures that 
such a mapping is invertible, i.e., that from xt the corre- 
sponding starting point xq can be uniquely determined, 
Xo = (/>-t(a;t). 

Consider now a general invertible and differentiable 
mapping 



(6) 



that maps phase space point x into phase space point x' . 
Here, the mapping (j){x) takes the place of the propagator 



(j)t{x). For such mappings Jarzynski has derived an ex- 
pression akin to the non-equilibrium work theorem |l2l |. 
To introduce the necessary notation we rcderivc this re- 
sult. For this purpose we consider the definition of the 
free energy difference, 

cxp{~PAF) = — = ^ . (7) 

WA Qa 

Multiplying and dividing the integrand in the above 
equation with exp{— /37i(^^^(x'), A^)} we obtain 



exp(-/3Ai^) 



Qa 



exp{-/3[H(x', As) - i(x'), A^)]}. (8) 

A change of integration variables from x' to x = (j>^^{x') 
yields 



exp(-/3AF) 



dx 



d(j){x) 



dx 



exp{-/3'H(x, Aa)} 
Qa 



cxp{~mWx), \b) ~ nix, \a)]}. (9) 

where \d(j){x)/dx\ is the Jacobian determinant of the 
mapping 0(x). The above equation suggests generaliz- 
ing the definition of work 

d(j){x) 



= H((/)(x), As) - H{x, Xa) - ksTln 



dx 



(10) 



This "work" includes the energy change caused by 
switching the control parameter from Aa to A^ . In addi- 
tion includes a term involving the Jacobian of 4'{x), 
which can be viewed as the work necessary to compress 
or expand the phase space volume when applying the 
mapping (/>(x). This entropic contribution can be inter- 
preted as "heat" absorbed during the mapping. In fact, 
if we choose the mapping to be the system's propaga- 
tor, this term is exactly the heat. In this interpretation, 
Equ. (|10|) is nothing other than an expression of the first 
law of thermodynamics. 

Using the work definition l|l()(l we can rewrite Equ. © 

as 



exp(-/3AF) ^ J dx 



exp{-/37i(x, Aa)} 
Qa 



!xp{-/3VU4x)} 
(11) 

or, as an average over the initial equilibrium distribution, 

(12) 



exp(-/3AF) - (exp{-/3W^(x)}). 

This equation can be viewed as a generalization of 
Jarzynski's identity. If the mapping is chosen to be the 
propagator (/)t{x) of Newtonian dynamics, the work 
equals the physical work W carried out on the system as 
it evolves from xq to xt, 

w = nixt,XB)-nixo,XA). (13) 

and Equ. (|12|) reduces to Jarzynski's identity. In deriv- 
ing this result we have exploited the fact that Newtonian 
dynamics conserves phase space volume - even when a 
control parameter changes with time, the Jacobian ap- 
pearing in the definition of [Equ. H1U|) ] is unity. 
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B. Long time step trajectories 



dn 



(16) 



Instead of the propagator (j>t , we can choose a sequence 
of molecular dynamics steps as our mapping. Each of 
these steps, which are designed to approximate the time 
evolution of the system over a small time interval At, 
maps a phase point Xi into a phase point Xi+i. Equa- 
tion (|12|l can be applied to a map defined by n such 
steps, together taking the initial point xq into a final 
point Xn = (f>n{xo)- The expression for the work is 
particularly simple for integrators such as the Verlct al- 
gorithm, which both conserve phase space volume and 
are time-reversible 0, E|- In this case the Jacobian of 
the mapping is unity, and, according to Equ. (|10|l . 



W^{xo) = H{xn, Xb) - Hixo, Aa), 



(14) 



so that 



-PAF 



(exp{-/3[7^(a;„. As) - 7^(xo, A^)]}) . (15) 



This relation is exact regardless of the size of the time 
step At used in applying these integrators. 

Equation (|15|l suggests the following algorithm. (1) A 
canonical distribution p{xo) of initial conditions is sam- 
pled with a Monte Carlo procedure or with an appro- 
priately thermostatted molecular dynamics simulation. 
(Note that in the latter case a sufficiently small time step 
must be used in order to preserve the correct equilibrium 
distribution.) (2) These initial conditions are then used 
as starting points for fast switching trajectories obtained 
by repeated application of the Verlet algorithm. (3) Dur- 
ing the integration the control parameter is changed from 
Xa to Xb- Since Equ. (|15|) is exact for any size of the 
time step At, the chosen integration time step can be 
arbitrarily large, provided that the variables specifying 
the state of the system (for instance, positions and mo- 
menta of all particles) retain values that do not exceed 
the range a computer can represent. We call this limit 
the stability limit. For each trajectory the energy dif- 
ference W = 'H{xn,XB) — 'H{xo,Xa) is determined and 
used to calculate the exponential average appearing in 
Equ. (|15|l . Since large time step trajectories are compu- 
tationally less expensive, this algorithm holds promise to 
increase the efficiency of fast-switching free energy calcu- 
lations. Whether this is actually the case depends on how 
the work distribution is modified by the increase in time 
step length. In Sec. Ill we describe how to analyze the 
efficiency of fast switching simulations with large time 
steps. 



C. Stochastic dynamics 

Often the dynamics of model molecular systems evolve 
by stochastic equations of motion. Common examples 
include the Langevin equation [T^ . 



where 7 is a friction coefficient and r]{t) is a fluctuat- 
ing random force; and deterministic dynamics coupled to 
stochastic thermostats, such as the Andersen thermostat 
[T^ . It has been shown that the Jarzinsky relation re- 
mains valid also in these cases 0, 0|. In this section 
we discuss the question whether the large time step ap- 
proach discussed in the previous section can be applied 
to stochastic dynamics as well. 

We describe the stochastic component of these dynam- 
ics through a "noise history" T]{t). In the case of Langevin 
dynamics the noise history is, as the notation suggests, 
the trajectory of the random force. For a given realiza- 
tion of r](t), the time evolution of a stochastic system can 
be regarded as deterministic, and we may write 



xt = <l>[xo;ii{t)]. 



(17) 



where the second argument in the deterministic map in- 
dicates its dependence on the noise history 77 (i). Since 
this mapping is invertible and differentiable, Equ. H12|) 
applies for any particular noise history, 

exp(-/?AF) = J dxpix)cxp{~(3W4x;r]it)]}), (18) 

where we have averaged over canonically distributed ini- 
tial conditions and 



W4x;v{t)] = 

H{<l>[x;Tj{t)],XB)~n{x,XA) - fcsTln 



dcl,[x;r,{t)] 



dx 



.(19) 



This result above is valid for any noise-dependent map 
4>[xq; ri{t)], so we are free to choose a mapping comprised 
of repeated application of Brownian dynamics steps 
with arbitrary step size. Because the result of averaging 
eyi'p{—(3W^[x;r](t)]} over initial conditions is completely 
independent of r](t), the remaining average over noise his- 
tories is trivial, yielding 

exp(-/3AF) = 

VT^it) f dxP[7^itMx)c^p{-pW4x;r,m-m 



Here, the notation J "Driit) indicates summation over all 
noise histories and P[i]{t)\ is the probability distribution 
for observing a particular realization. 

Interestingly, the above derivation implies that Equ. 
(|20|) can be applied to mappings with a completely ar- 
bitrary stochastic component. In particular, it is not 
necessary that the magnitude of stochastic fluctuations 
be related in any way to the rate of dissipation. In order 
to preserve a canonical distribution, the Langevin equa- 
tion must be supplemented with such a constraint on 
the statistics of rj{t) (henceforth assumed to be Gaussian 
white noise): 



q = p/j 



{■q{t)'n{t'))=2kjiT-f5{t~t'). 



(21) 
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This fluctuation-dissipation relation, ensuring detailed 
balance, is a necessary condition for the applicability of 
Jarzynski's original identity concerning the exponential 
average of work defined in the conventional way. Obtain- 
ing that result from Eau. I20| is not nearly as straightfor- 
ward as was the analogous task for deterministic dynam- 
ics. For stochastic mappings the Jacobian determinant 
does not directly correspond to heat, even in the limit of 
small time step size. Instead, for a trajectory of length r 



lim 



dx 



(22) 



where rif is the number of momentum degrees of free- 
dom. The volume of a phase space element evolving 
under Langevin dynamics thus decays steadily as time 
evolves and has no contribution from the fluctuating ran- 
dom force. Mathematically, this phase space compres- 
sion arises from the systematic damping of kinetic en- 
ergy through the friction term in Equ. El Physically, 
its cancellation of contributions from heat Q in the ex- 
ponential average of Equ. 1201 is a subtle consequence of 
detailed balance. Crooks has shown that, for trajecto- 
ries generated by any balanced dynamical rules, e~^'^ is 
equivalent to the ratio of probability densities of forward 
and time-reversed pathways. This ratio is closely related 
to a mapping's Jacobian, as will be shown in Sec. |V] 

The practical utility of our large time step result for 
free energy differences is compromised in the specific case 
of Langevin dynamics by at least two issues. First, the 
rapid decay of \d4>[x\rj{t)]/ dx\ could damp all but the 
largest fluctuations, making convergence of the average 
in Equ. I2UI problematic. Second, an exact calculation 
of the Jacobian for large time steps is cumbersome when 
many degrees of freedom interact. The insensitivity Equ. 
(|20|1 to the form of stochastic noise might be exploited to 
offset these problems, but it is not obvious how to do so. 



III. EFFICIENCY ANALYSIS 

In a straightforward application of the fast switching 
procedure, the free energy difference Ai^ is estimated 
from a finite sample of N trajectories originating from 
canonically distributed initial conditions. Here we as- 
sume that these initial conditions are statistically uncor- 
related samples. Defining 

X = exp(-/?iy^), (23) 

we can write the free energy difference for finite TV as 

AFat EE -/cbTIuXa,. (24) 

where Xat is the average of X over N independent tra- 
jectories: 



1 ^ 



N 



(25) 



Here, X^*) is the value of X obtained from the i-th tra- 
jectory. We now repeat this entire procedure M times 
(generating a total of MN trajectories) and average over 

all M resulting free energy estimates Ai^^ . The limiting 
result of this protocol is 



N 



M 

hm ^VAFj^'. 



(26) 



Due to the non-linearity of the logarithm relating the 
average to the free energy estimate (see Equ. 
(|24|l '). (AFff) differs from the true free energy difference 
AF even in the limit of infinitely many repetitions. For 
sufficiently large N this deviation, or bias, is given by 



5Ar 



{AF 



N 



AF 



knT {{5Xf 
2N (X)2 



(27) 



where the fluctuation 5X = X — {X) is the deviation 
of X from its average value. This equation is obtained 
by expanding the logarithm in Equ. (|24|l in powers of 
relative fluctuations of Xjv and truncating the expansion 
after the quadratic term. 

In addition to the systematic bias, we must account 

for random error in each free energy estimate F^^ . As- 
suming the statistical errors of different samples to be 
uncorrelated, and denoting their variance as 



nl = {[AFn - {AF 



N 



(28) 



we obtain the total mean squared deviation from the true 
free energy difference. 



2 



\AF 



N 



AFl 



2 



(29) 



Note that the bias decays with sample size as 6jv oc 1 /N 
for large N, while the scale of random error decays only 
as (Tat cx 1/^/N. Thus, for sufficiently large N only the 
statistical errors in the free energy estimate are relevant, 
and we can safely approximate 



fc|T2 



{{SX)' 



'N 



N {xy 



(30) 



From Equ. (EOJ 

one can determine the number of tra- 
jectories necessary to obtain a certain level of error 



{X? 



(31) 



The computational cost of each trajectory is roughly pro- 
portional to the number of required force calculations and 
hence proportional to the number of steps n = t j At nec- 
essary to generate a trajectory of length r. Neglecting 
the cost of the generation of initial conditions we can thus 
define a normalized computational cost 



Ccpu = 



_{m?) 



1=1 



{xy 



At {xy 



(32) 



5 



This computational cost Ccpu is the CPU-time required 
to obtain an accuracy in the free energy of e = k-gT 
measured in units of the CPU-time required to carry out 
one single molecular dynamics time step. Note that while 
{X) is independent from the stcpsize At, the mean square 
fluctuations {{SX)'^) depend on it. To determine an op- 
timal time step size for a fast-switching free energy cal- 
culation, we must minimize the entire quantity Ccpu- In 
the following section we will present calculations of this 
normalized CPU-time as a function of the stepsize At for 
two different models. 



IV. NUMERICAL RESULTS 

A. One-dimensional system 

We first study the effect of large time step integration 
for a simple model introduced by Sun 0- In this one- 
dimensional model a point particle of unit mass moves in 
a potential that depends on a control parameter A. The 
Hamiltonian for the system as a function of the position 
q and the momentum p of the moving particle is 



n{q,p; X) ^ -p^ + q""- 16(1 -X)q^, 



(33) 



Here, all quantities are scaled to be unitlcss. For A = 
the potential has two symmetric minima located at 
separated by a barrier of height AE = 64 
located at q = 0. For A = 1 the potential energy 
function reduces to a single quartic well. For a given 
value of the control parameter A the partition function 
is / dqdpexp{—PT-l{q,p,X)}. The free energy difference 
between the two states corresponding to A = 1 and 
A = 0, respectively, can be calculated analytically to be 
62.9407fcBr [ni. 

For this model we carried out fast switching simula- 
tions using different time steps At ranging from At = 
0.002 to At = 0.1. The equations of motion were inte- 
grated with the velocity Verlet algorithm ^| , yield- 
ing positions qt and momenta pt as a function of time 
t. For this model the Verlet algorithm becomes unsta- 
ble for time steps larger than At = 0.1. In all cases the 
total trajectory length was t = 10, corresponding to a 
transformation sufficiently gradual to allow accurate cal- 
culation of the free energy difference and of the mean 
square fluctuations {{SX)'^). For each time step size we 
integrated 10^ trajectories from initial conditions sam- 
pled from a canonical distribution (with A = 0) using a 
Monte Carlo procedure. Along the trajectories the con- 
trol parameter was varied from its initial to its final value. 
More precisely, after each velocity Verlet step, carried 
out at constant A, the control parameter was increased 
by l/n where n = r/At is the number of time steps in 
the trajectory. At the end of each trajectory the work 
W = H{qr,PT'A) — 'H{qo,Po',0) was calculated and added 
to the exponential average. (Note that, when advancing 
time in large steps, it is important to use this generalized 



definition of work rather than summing estimates of the 
physical work performed during each stepwise change of 
the control parameter.) 




FIG. 1: Work distributions P{W) obtained for the Sun model 
for a trajectory length t — 10 and step sizes At = 0.1 (solid 
line) and 0.01 (dashed line). Work distributions for all time 
steps smaller than At = 0.01 are indistinguishable from the 
distribution for At — 0.01 on the scale of the figure. 

Work distributions P{W) obtained for different step 
sizes At (and hence for different numbers of steps per 
trajectory) are shown in Fig. ^ P{W) deviates visibly 
from its small time step limit only for the largest step 
size. At = 0.1 These differences originate in the inac- 
curacy of the integration algorithm for large step sizes. 
Even though the work distribution varies with the size 
of integration steps, the resulting free energies show no 
step size dependence (see Fig. E)), provided the stability 
limit of the integration algorithm is not exceeded. 
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FIG. 2: Free energy differences (circles) obtained for the Sun 
model from fast switching trajectories of length r = 10 with 
different step sizes At. The dotted line denotes the exact 
free energy difference. Also shown is the average work (W) 
(squares). 

To quantify the statistical error in the free energy es- 
timates shown in Fig. [2 we have calculated the relative 
fluctuations {{5 X^) / {X)"^ , which according to Equ. l|Sn|) 
determine the mean squared error e^. These relative 
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fluctuations are plotted in Fig. |21 The irregular shape 
of this curve reflects changes in the features (such as the 
peak near W = 62) of the work distributions shown in 

Fig. HI 
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FIG. 3: Relative fluctuations {{SXf)/{Xf for the Sun model 
as a function of the step size At. 



The computational cost Ccpu follows from the relative 
fluctuations and is shown as a function of the step size 
in Fig. ^ Over the range of time steps depicted in the 
figure the computational cost decreases from about 10^ 
to about 10^. Thus, the fast switching simulation can 
be accelerated by two orders of magnitude if the conser- 
vative step size of At = 0.001 is replaced by a step size 
At = 0.1 near the stability limit. Although in the lat- 
ter case the equations of motion are not faithfully solved 
(see Sec. Fig. I10|l . the expression for the free energy 
remains exact. 
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FIG. 4: Normalized CPU time Ccpu for the Sun model as a 
function of the step size At for ksT — 1. For all step sizes 
the total trajectory length was r = 1. This result indicates 
that for this model the computational cost of a free energy 
calculation with a given error decreases for increasing step 
size until the stability limit is reached. 



B. Dragged particle in Lennard- Jones fluid 

Our second example involves a system with many de- 
grees of freedom and is therefore likely to be more rele- 
vant for typical molecular systems of interest. Spepcifi- 
cally, we tested the large-timestep version of the Jarzyn- 
ski identity for a particle dragged through a Lennard- 
Joncs fluid. In these simulations, a tagged particle is 
coupled to a harmonic trap whose minimum is shifted 
from one position to another while the system evolves in 
time. If this process is carried out at a finite rate, work 
is performed on the system by the moving trap. Never- 
theless, the free energy difference between the two states 
corresponding to the initial and final position of the trap 
vanishes due to symmetry. 

The time-dependent Hamiltonian for our Af-particle 
system is 

M 2 I. 

H{p, 9, = E ^ + E "(^^^ ) + o - (34) 



i=l 



where Pi is the momentum of particle i, is the position 
of particle i, Tij is the distance between particles i and j 
and the second sum on the right hand side extends over 
all particle pairs. The particles interact via the Lennard- 
Jones potential 



v[r) = 4e 



(35) 



Here, e and a are parameters describing the depth of the 
potential well and the interaction range of the strongly 
repulsive core, respectively. The last term in Equ. Htj4|) 
describes the potential energy of one particular parti- 
cle (we arbitrarily pick particle 1) in a harmonic trap 
with force constant k. In Equ. R(t) denotes the 

time-dependent position of the trap's minimum, which is 
moved from the origin in the x-direction with constant 
speed V, 



R(t) = e^ji/t, 



(36) 



where e^; is the unit vector in cc-direction. During a to- 
tal time r the trap is displaced by an amount L = vt. 
While the trap's minimum roughly determines the posi- 
tion of particle 1. this particle does fluctuate about R(t). 
The work W performed on the system along a particular 
trajectory from {(7o,Po} to {qTjPr} is given by 



= H{qr,PT\r) - 7i{qa,po;0)- 



(37) 



Since the free energy of the system does not depend on 
the trap's location, AF = 0, the exponential work av- 
erage carried out over many realizations of this process 
is 



(exp(-/3W^)) ^ cxp{-PAF) = 1. 



(38) 



Here, angular brackets indicate a canonical average for a 
fixed position of the trap. 
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We have carried out a fast switching procedure for 
M = 108 particles of unit mass in a three-dimensional, 
cubic simulation box with periodic boundary conditions. 
The fluid's density, p — 0.8cr~^, is roughly that at 
the triple point, and its initial temperature places it in 
the liquid phase. Canonically distributed initial condi- 
tions were generated by a molecular dynamics simula- 
tion, employing an Andersen thermostat at temperature 
k-gT/e = 1.0 [iJI, with the trap fixed at the origin. For 
the generation of initial conditions, the equations of mo- 
tions were integrated with the velocity Verlet algorithm 
and a time step of dt = 0.001cr(e/m)-'^/^. The state of the 
system was recorded every 50 steps along this equilibrium 
trajectory, providing an ensemble of intial conditions for 
the fast switching procedure. 

From these initial conditions we generated fast switch- 
ing trajectories of total length r = 1.2a{e/m)^/^ with 
different time steps ranging from At = 0.00lCT(£/m)^/^ 
to At = 0.02(T(e/m)^/^. The corresponding pathways 
ranged in number of steps from n = 1200 to n = 60, 
respectively. In each case the total displacement of the 
trap from its initial to its final position was L = 0.5cr 
corresponding to a velocity of v — ^jYlimj e)^!'^ . The 
velocity Verlet algorithm without thermostat was used 
to integrate Newton's equations of motion. Along these 
trajectories the particle trap was displaced stepwise by a 
small distance Ljn after each molecular dynamics step. 
Work distributions obtained in this manner are shown in 
Fig. |S1 for three different step sizes. Free energy differ- 
ences AF and the average work W calculated in these 
simulations are depicted in Fig. El 
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FIG. 5: Work distributions PiyV) for a particle dragged 
through a Lennard- Jones fluid for step size At = 
0.001 cr(e/m)^''^ (solid line), step size At = 0.015 cr(e/m)^''^ 
(dashed line) and step size At = 0.02(T(£/m)^/^ (dotted line). 

We used the relative fluctuations {{5XY) / {X)"^ (see 
Fig. O to calculate the normalized CPU time Ccpui 
which is plotted in Fig. |S1 The computational ef- 
fort required to obtain a specific accuracy decreases 

with increasing step size until the largest value. At ~ 
0.02cr(£/m)^/2_ jygt 

as in the schematic one-dimensional 
example, the most efficient fast-switching calculation for 
dragging a particle through a dense fluid is obtained for 



< 



2 - 



□ u- 



Oh o- 

I III 



-o- 

I 



J 



0.001 



0.01 



At 



FIG. 6: Free energy difl'erences (circles) for a particle dragged 
through a Lennard- Jones fluid as a function of step size At. 
The trajectory length was r = 1.2(j{e/mY^'^ for all step sizes. 
Also shown is the average work {W) (squares). 



step sizes close to the stability limit. 
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FIG. 7: Relative fluctuations {{5Xf)/{Xf for the particle 
dragged through a Lennard- Jones fluid as a function of the 
step size At. 



V. DISCUSSION 

Molecular dynamics simulations carried out with large 
time steps do not faithfully reproduce the dynamics of 
any system with configuration-dependent forces. If such 
large time steps are used in the generation of fast switch- 
ing trajectories to calculate free energy differences on the 
basis of Jarzynski's identity, integration errors lead to 
work distributions differing from those obtained in the 
small time step limit. Nevertheless, as we have shown 
in Sec. ^1 the Jarzynski identity remains exactly valid 
in principle for time steps of arbitrary size. As a prac- 
tical matter the stability limit provides an upper bound 
to the step size. Since the computational cost of molec- 
ular dynamics trajectories is proportional to the number 
of integration steps, large time steps can be beneficial. 
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FIG. 8: Normalized CPU time Ccpu for the particle dragged 
through a Lennard- Jones fluid of M = 108 particles at a den- 
sity p = 0.8a~^ as a function of the step size At at ksT/e = 1. 
The particle is dragged through the fluid by a parabolic po- 
tential with force constant k — 10''(e/a^) moving at constant 
speed u — 5/12(m/e)^'''^. For all step sizes the total trajec- 
tory length was t = 1.2(T(e/m)^''^. Also for this model the 
cost of a free energy calculation decreases up to step sizes just 
short of the stability limit. 



Whether an increase of computational efhcieny can re- 
ally be achieved depends on how the relative squared 
fluctuations {{SX)'^) / (X)'^ scale with the size of the inte- 
gration time step. If these fluctuations increase with At 
sublinearly, using large time steps is advantageous. 

In both of the numerical examples presented in Sec. lIVI 
increasing the size of the time step well beyond the range 
appropriate for equilibrium simulations proved favorable. 
Resulting efficiency increases reached up to two orders 
of magnitude. Here we ask more generally, and with 
complex molecular systems in mind, what circumstances 
should allow large time step integration to improve upon 
fast switching simulations. To answer that question it is 
convenient to write the work W performed on the system 
during the transformation as sum of two parts. This 
separation is particularly natural if the control parameter 
A is changed stepwise, i.e., if each step is performed at 
constant A and is followed by an increase in the control 
parameter by an increment AA = 

We define the first contribution to the generalized work 
as the change in energy due to the changes in control 
parameter for fixed phase points. 



^7i[a-,,iAA] 1)AA]. (39) 



i=l 



Here, Xi is the phase space point reached after i integra- 
tion steps starting from phase space point xq. In Crooks 's 
considerations of systems out of equilibrium, Wx is pre- 
cisely the physical work exerted by the change of control 
parameter 

In a system that evolves according to Newton's equa- 
tions with a time-independent potential energy function, 
the total energy is a constant of the motion. But when 



these equations are integrated approximately over finite 
time steps, the energy of the system is not perfectly con- 
served. Summing the energy changes due to integration 
error in the intervals between stepwise increases of A, we 
obtain the other contribution to W^: 



E 

i=0 



H[xi+i,iA\] — T-l[xi, iAX]. 



(40) 



To summarize our decomposition of W^, the work 
Wx involves changes in control parameter at fixed phase 
space points, while the less physical contribution We in- 
volves changes in the phase space point at constant con- 
trol parameter. The total work performed during the 
transformation is the sum of these two quantities, 



W = Wx + We^ 



(41) 



In the limit of very small time steps, the energy of the sys- 
tem is conserved whenever A is constant. In this case, the 
"error" work W^ vanishes and the entire work is caused by 
changes in the control parameter, W = Wx- If a system 
is driven away from equilibrium by a stepwise increase 
of the control parameter, the error work serves as a sim- 
ple measure for the accuracy of approximate numerical 
integration. 

Neglecting correlations between Wx and W^ for the 
moment, we can write 

(exp(-/3AVF)) - (exp(-/31^A))(cxp(-/3W,)). (42) 

Accordingly, the free energy change is the sum of two 
terms originating from Wx and We, 



AF = AFx + AFe 

= -fcBrin(e-''^^ 



fcBnn(e 



(43) 



This separation of the free energy difference into two 
terms related to Wx and We , respectively, is only strictly 
valid if fluctuations in 



Xx = eyi-p{~f3Wx) and Xe = eyi-p{-(3We, 



(44) 



are statistically independent. This supposition cannot 
be assumed a priori to be the case. For the two models 
treated numerically in this study we have calculated the 
correlation 



Cxe = 



{SXxSXe 



(45) 



where SXx = Xx - (Xx) and SXe = Xe - {Xe). This 
coefficient quantifies correlations between the two vari- 
ables Xx and Xe- While in the absence of correlations 
Cxe = 0, perfect correlation (anticorrelation) leads to 
Cxe — 1 {Cxe = ~1)- Correlations computed numerically 
for the Sun model are depicted in Fig. |51 as a function of 
the time step A^. For all time step sizes Xx and Xe are 
only weakly (anti)correlated. The assumption of statis- 
tical independence is justified in these cases. 
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FIG. 9: Correlations C^x between control parameter work and 
integration error work as a function of the step size for the 
Sun model. 



Under the same assumption the statistical error of a 
fast switching free energy calculation with large time 
steps can also be written as the sum of two distinct contri- 
butions. Absent correlations between W\ and W^, Equ. 
becomes 



N 



2 r 



(46) 



Thus, the total mean squared error is the sum of the 
mean squared errors of the free energies related to the 
control parameter work and the integration error work, 
respectively. One potentially substantial difference be- 
tween these two error contributions is their dependence 
on system size. Often, the control parameter A acts only 
on a small subset of the system's degrees of freedom. In 
calculating the chemical potential of an electrically neu- 
tral species by particle insertion, for instance, the solute 
interacts only with a small number of other particles near 
the insertion point. Similarly, the transformation of a 
residue of a protein into another one mainly affects only 
a local group of interactions. In such cases, the work W\ 
resulting directly from changes in the control parameter 
quickly saturates with growing system size. The work 
Wf: , originating in the inaccuracy of the integration algo- 
rithm, has a different system size dependence. Since all 
degrees of freedom contribute to the integration error, We 
is expected to grow linearly with the system size. Thus, 
for sufficiently large systems the second term on the right 
hand side of Equ. H4t)|) might become dominating for 
long time steps, thus limiting the maximally possible ef- 
ficiency gain. (It is for similar reasons that hybrid Monte 
Carlo-molecular dyiiamics simulations decline efficiency 
for large systems But as long as the integration er- 

ror is small compared to work done by changes in control 
parameter, using large time steps should remain advan- 
tageous. 

In our discussion of the efficiency of the large time step 
fast switching approach we have so far neglected the com- 
putational cost associated with generating initial condi- 
tions. Often, starting points for the fast switching tra- 



jectories are generated with an appropriately thermostat- 
ted molecular dynamics simulation. It is important that 
these simulations are carried out with a time step of con- 
ventionally small size. Otherwise the distribution of ini- 
tial conditions can differ from the necessary canonical 
one. Thus, the large time step approach can reduce only 
the computational cost associated with generating the 
non-equilibrium trajectories. This latter contribution is 
dominant by far in most cases. 

The issues discussed in this paper have some inter- 
esting implications for molecular dynamics simulations 
carried out at equilibrium as well. If the control param- 
eter A is not changed as the system evolves in time, no 
work W\ is done on the system. In this case the free 
energy difference AF also vanishes, since initial and final 
Hamiltonians are identical. Nevertheless, due to the im- 
perfection of a finite time step integration algorithm, the 
energy is not strictly conserved, and some "error" work 
We is performed. According to Equ. (|15|l the statistics 
of such energy errors obeys 



(eM-PWe)) = 1. 



(47) 



This result implies that the distribution of integration 
errors cannot be symmetric about W^ = 0. Rather, pos- 
itive errors are in loose terms more likely than negative 
ones. For a Gaussian error distribution, Equ. (|47|l relates 
the average error to the width of the error distribution, 



1 



(We) = -P'a, 



(48) 



For more complicated error distributions, this same re- 
sult is obtained from a cumulant expansion of Equ. H47|) 
truncated after the second term. Thus, for an error distri- 
bution that is Gaussian and/or sufficiently narrow, the 
average error is positive. This demonstration that the 
average energy of a simulated canonical ensemble drifts 
upward with time makes no reference to the details of 
molecular interactions. 

We have verified Equ. (|47|l both for the Sun model 
and for the Lennard-Joncs fluid. In both cases we have 
determined the distribution P{We) of the integration er- 
ror We for a large integration step At and constant A. 
The averages were carried out in the respective canonical 
ensembles corresponding to A = 0. Two typical distribu- 
tions are shown in Fig. ^| For the Sun model (solid line) 
P{We) is strongly non-Gaussian and asymmetric about 
We = 0, with positive errors more likely than negative 
ones. While the average error of (We) = 0.041 is clearly 
positive, the exponential average (exp(— /3We)) is unity 
with high accuracy. For the Lennard- Jones fluid with 
one particle in a fixed trap the distribution of integra- 
tion errors is approximately Gaussian, with an average 
of (We) ~ 0.064. Also in this case the exponential work 
average is unity. 

Formulating changes in control parameter as a se- 
quence of discrete steps is also convenient for demonstrat- 
ing the equivalence between our identity for stochastic 
mappings, Equ. (|20|1 . and Jarzynski's original identity. 
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FIG. 10: Distributions P{Wt) of the integration error Wc 
at constant control parameter. Solid line: Sun model for 
A = 0, r = 10, and At = 4/30; dashed line: Lennard- 
Jones fluid with one particle in fixed trap, r = 1.2(7(e/m)^'^^, 
At = 0.015 cr(e/m)^/^ 



During periods when A is fixed, stochastic evolution sub- 
ject to a suitable fluctuation-dissipation relation, such as 
Equ. H21(l . satisfies detailed balance in the small time 
step limit. Specifically, 



Pi{x)p{x — > x) = pi{x')p{x' — > x) 



(49) 



where pi{x) oc exp[— /37i(a;, iAA)] is the canonical distri- 
bution corresponding to the Hamiltonian Ti.{x,iAX) at 
step i, and pi{x ^ x') is the noise-averaged probability 
for a phase space point x at the beginning of a constant-A 
interval to evolve into phase space point x' at the end of 
the interval. The caret in Equ. H49|) indicates time rever- 
sal, so that p{x' x) is the probability for x' to evolve 
into X under dynamics running backward in time. For 
Langevin dynamics, time reversal can be achieved sim- 
ply by inverting the signs of all momenta contained in 
phase space points x and x' . For a given noise history, 
such transition probabilities are specified by the deter- 
ministic map (j)[x;ri{t)]. We can thus rewrite Equ. 



-pni.^AX) (^S{x'-^[x;^m)^ 



dcl)[x;v{t)] 



dx 



,(50) 



where angled brackets with a subscript r] denote an av- 
erage over realizations of the noise history. 

From the condition of detailed balance, we now obtain 
a general identity for averages involving the Jacobian de- 
terminant. We begin by multiplying Equ. (|50|) by an ar- 
bitrary function f(x,x') of two phase space points, and 
then integrate over x', yielding 



(/(x,0[x;,7(i)])), = (/(x,0[x;ry(t)]) 



xe 



-i3ini<pixi7]{t)],iA\)-nix,iA\)] 



dx 



(51) 



This result holds for any period during which the con- 
trol parameter is held constant, i.e., while no physical 
work is done on a system. Since we have assumed here 
that equations of motion are integrated with an infinites- 
imally small time step, the "error" work vanishes as 
weU. The energy change H{(l)[x] ry(t)], iAA) — Tl{x, iAX)] 
in Equ. (|?T)l is thus identically the heat Qi absorbed from 
the bath during the zth time interval. Equation l|51|l 
states that the Jacobian determinant and e~^^ negate 
one another in averages over noise history. 

This proof is completed by decomposing the total 
change in energy along a trajectory into contributions 
from physical work W\ (accumulated over many discrete 
steps in control parameter at fixed phase space point) 
and heat (accumulated over many intervals during which 
the phase space point evolves at fixed A): 

niq^[x; vit)lXB) - nix, Xa) ^Wx + Q (52) 

n-l 

Q = Y,n[x,+i,iAX]-H[x„iAX] (53) 

i=0 

We have retained the definition of physical work from 
Equ. (|40|l . Note that heat in the case of stochastic dy- 
namics is defined in precisely the same way as "error" 
work was defined for deterministic dynamics propagated 
in an approximate way. Recalling the generalized defi- 
nition of work and applying the identity (|51|l with 
f{x,x') = exp{-l3{H[x',iAX] - n[x',{i - 1)AA])} for 
each constant-A interval, we finally have 



d(f>[x;7]{t)] 



dx 



(54) 



This equivalence, together with Equ. (|20() . completes 
an alternative route to Jarzynski's identity for systems 
evolving stochastically under the constraint of detailed 
balance. 

The consequence of detailed balance expressed in this 
way has interesting implications for energy fluctuations of 
stochastic systems at equilibrium (i.e., with fixed control 
parameter). With the choice f{x,x') = 1, we have for 
the specific case of Langevin dynamics (see Equ. (|221l) 



(55) 



along trajectories of length r. This identity stands in 
stark constrast to the corresponding exponential average 
of energy fiuctuations under norm-conserving determin- 
istic mappings, Equ. H47|l . The long-time divergence in 
Equ. H55|) would be expected for a system which asymp- 
totically loses all memory of its initial conditions. In that 
case, the exponential average factorizes in the long-time 
limit: 



(56) 



The first factor on the right hand side of Equ. H56|l av- 
erages a quantity that negates the effect of Boltzmann 
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weighting, and is proportional to the entire phase space 
volume. Since the range of possible momenta is unbound 
even in a system with finite volume, a divergent result is 
inevitable once correlations have decayed completely. We 
anticipate similarly unbounded growth of exponentially 
averaged energy fluctuations for many classes of stochas- 
tic dynamics, such as Monte Carlo sampling. That the 
analogous average is fixed at unity for deterministic prop- 
agation rules with unit Jacobian such as the Verlet algo- 
rithm indicates that errors arising from finite time step 
size do not disrupt substantial correlations with initial 
conditions. 



VI. CONCLUSION 

By considering general invertible phase space map- 
pings we have demonstrated that the Jarzynski rela- 
tion remains exactly valid for non-equilibirum trajecto- 
ries generated with large time steps, provided the work 
performed on the system is defined appropriately. For in- 
tegration algorithms that conserve phase space volume, 
such as the Verlet algorithm, this definition is partic- 
ularly simple. Here, the work just equals the energy 
difference between the final and the initial state of the 
trajectory. Simulating dynamics with a larger time step 



requires fewer integration steps to generate a trajectory 
of given length, and therefore lower computational cost. 
Numerical simulations indicate that optimum efficiency 
is achieved for time steps just short of the stability limit. 
Compared to simulations with time steps of conventional 
size, the long time step approach can yield improvements 
in efficiency of one or more orders of magnitude. 

Recently, Sun has shown how work-biased path sam- 
pling can be used to improve the efficiency of fast switch- 
ing simulations 0. However, it seems that this path sam- 
pling approach does not outperform conventional meth- 
ods for calculating free energy differences. It will be in- 
teresting to see if the fast switching approach can be im- 
proved, by combining the long time step approach of this 
paper with biased path samping methodologies |^l3.llll|. 
to the point that it is computationally competitive with 
other free energy calculation techniques, such as umbrella 
sampling, thermodynamic integration, or flat histogram 
sampling. 
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